import os
from pyx import *
from pyx.deco import barrow,earrow
from pyx.style import linewidth, linestyle
from pyx.graph.axis import painter, tick
from pyx.graph.axis import *
from scipy import *
from scipy.optimize import leastsq

unit.set(xscale=2)

def peval(x, p):
	return p[0]-p[1]*x**p[2];
#       return p[0]-p[1]*x;

def fit(x,y):
#    p0 = [1, 1,2]
    p0 = [1, 1, 1]
    plsq = leastsq(residuals, p0, args=(y, x))
    print plsq[0]
    return x,plsq[0]

def residuals(p, y, x):
    p1,p2,p3 = p
#    err = y-p1+p2*x**p3
    err = y-p1+p2*x
    return err


os.system('rm dat.txt')
i=4
while i<100 :
	os.system(('./el3 --N '+str(i))) 
	i=i*2




data=io.array_import.read_array('dat.txt')
NbDeTest=len(data[:,0])
xf1,yf1 = fit( log(data[:,0]),log(data[:,1]))
print xf1,yf1,NbDeTest
d=[[x,peval(x,yf1)] for x in xf1]
slope1=yf1[1]


xf2,yf2 = fit( log(data[:,0]),log(data[:,2]))
print xf2,yf2,NbDeTest
d=[[x,peval(x,yf2)] for x in xf2]
slope2=yf2[1]


g = graph.graphxy(width=29,x=graph.axis.log(title="$h$"),
                  y=graph.axis.log(),key=graph.key.key())


g.plot(graph.data.file("dat.txt", x=1, y=2,title="Norme H1"),styles=[graph.style.line([color.rgb.red])])
g.plot(graph.data.file("dat.txt", x=1, y=3,title="Norme L2"),[graph.style.line([color.rgb.green])])
g.text(g.width/2, g.height + 0.2, "Estimation de l'erreur en base loglog", 
       [text.halign.center, text.valign.bottom, text.size.Large])
g.finish()

x1, y1 = g.pos((data[NbDeTest/2+1,0]),(data[NbDeTest/2+1,1]))
x2, y2 = g.pos((data[NbDeTest/2-1,0]),(data[NbDeTest/2+1,1]))
x3, y3 = g.pos((data[NbDeTest/2-1,0]),(data[NbDeTest/2-1,1]))
g.stroke(path.line(x1-.5, y1, x2+.5, y2), [linestyle.dashed])
g.stroke(path.line(x2, y2-.5, x3, y3+.5), [linestyle.dashed])
g.stroke(path.line(x1, y1-.5, x1, y1+.5), [linestyle.dashed])
g.stroke(path.line(x3-.5, y3, x3+.5, y3), [linestyle.dashed])
g.stroke(path.line(x1,y1-.25,x2,y2-.25 ), [barrow.normal, earrow.normal])
g.stroke(path.line(x2+.25,y2,x3+.25,y3 ), [barrow.normal, earrow.normal])
g.text((x1+x2)/2,y1-.5, r"$a1$", [text.vshift.middlezero])
g.text(x2+.5,(y2+y3)/2, r"$b1$", [text.vshift.middlezero])
g.text(x1, y3, r"$b1/a1 = $"+str(abs(slope1)), [text.halign.center])

x1, y1 = g.pos((data[NbDeTest/2+1,0]),(data[NbDeTest/2+1,2]))
x2, y2 = g.pos((data[NbDeTest/2-1,0]),(data[NbDeTest/2+1,2]))
x3, y3 = g.pos((data[NbDeTest/2-1,0]),(data[NbDeTest/2-1,2]))
g.stroke(path.line(x1-.5, y1, x2+.5, y2), [linestyle.dashed])
g.stroke(path.line(x2, y2-.5, x3, y3+.5), [linestyle.dashed])
g.stroke(path.line(x1, y1-.5, x1, y1+.5), [linestyle.dashed])
g.stroke(path.line(x3-.5, y3, x3+.5, y3), [linestyle.dashed])
g.stroke(path.line(x1,y1-.25,x2,y2-.25 ), [barrow.normal, earrow.normal])
g.stroke(path.line(x2+.25,y2,x3+.25,y3 ), [barrow.normal, earrow.normal])
g.text((x1+x2)/2,y1-.5, r"$a2$", [text.vshift.middlezero])
g.text(x2+.5,(y2+y3)/2, r"$b2$", [text.vshift.middlezero])
g.text(x1, y3, r"$b2/a2 = $"+str(abs(slope2)), [text.halign.center])

g.writeEPSfile("graph")

